#####Measurement Error and the Specification of the Weights Matrix in Spatial
#####Regression Models by Garrett N. Vande Kamp

#####Replication File

#2. Clear the workspace and set a working directory
rm(list = ls())
setwd()

#3. Load packages
library(sem)
library(AER)
library(xtable)
library(multcomp)
library(Cairo)
library(extrafontdb)
library(extrafont)



#3.5 Time the Program
Start.Time <- Sys.time()

#4. Set the seed
set.seed(8675309)

#5. Define obKects for running the simulation
spectral <- function(W){
  Decomp <- eigen(W)
  W.normalize <- W/max(Decomp$values)
}

knn <- function(D, k, h){
  diag(D) <- h
  for (Z in 1:k){
    for (i in 1:n){
      Col.Min <- min(D[i,])
      D[i,] <-ifelse(D[i,]==Col.Min, h+1, D[i,])
    }}
  D<- ifelse(D==h+1,1,0)}

coverage <- function(b, se, true, level=.95, df=Inf){
  qtile <- level + (1-level)/2
  lower.bound <- b - qt(qtile, df=df)*se
  upper.bound <- b + qt(qtile, df=df)*se
  true.in.ci <- ifelse(true >= lower.bound & true <= upper.bound, 1, 0)
  cp <- mean(true.in.ci)
  mc.lower.bound <- cp - 1.96*sqrt((cp*(1-cp))/length(b))
  mc.upper.bound <- cp + 1.96*sqrt((cp*(1-cp))/length(b))
  return(list(coverage.probability=cp,
              true.in.ci = true.in.ci,
              ci = cbind(lower.bound, upper.bound),
              mc.eb = c(mc.lower.bound, mc.upper.bound)))
  
}
#5.5: Define Distance Matrix and Generate W Matrices 
reps <- 1000
obs <- 1000
n <- 100
t <- 10
I <- matrix(0, nrow=obs, ncol = obs)
diag(I) <- 1
x.coord <- runif(obs,0,100)
y.coord <- runif(obs,0,100)
points <- cbind(x.coord,y.coord)
D<- as.matrix(dist(points))
Wlinraw <- (max(D)-D)/max(D)
diag(Wlinraw) <- 0
Winvraw<-1/D
diag(Winvraw) <- 0
Winvsquaredraw <- 1/(D^2)
diag(Winvsquaredraw) <- 0
K.raw<- matrix(1, nrow=obs, ncol=obs)
diag(K.raw) <- 0


x.coord.small <- runif(n,0,100)
y.coord.small <- runif(n,0,100)
points.small <- cbind(x.coord.small,y.coord.small)
D.small<- as.matrix(dist(points.small))
I.T <- matrix(0, nrow=t, ncol=t)
diag(I.T) <- 1
D.panel <- kronecker(I.T,D.small)

matrixofones <- matrix(1, nrow=n, ncol=n)
Ones.T <- kronecker(I.T,matrixofones)

K.panel <- Ones.T
diag(K.panel) <- 0
Wlinraw.panel <- (max(D.panel)-D.panel)/max(D.panel)
Wlinraw.panel <- Wlinraw.panel*K.panel
Winvraw.panel<-1/D.panel
for(i in 1:obs) {
  for(K in 1:obs){
    Winvraw.panel[i,K] <- ifelse(Winvraw.panel[i,K]==Inf, 0, Winvraw.panel[i,K])
  }
}


#6: Define the systematic component of the model
X <- runif(obs, -1, 1)
Winvspec <- spectral(Winvraw)
Winvspec.panel <- spectral(Winvraw.panel)
Wlinspec <- spectral(Wlinraw)
Wlinspec.panel <- spectral(Wlinraw.panel)
Kspec <- spectral(K.raw)
Kspec.panel <- spectral(K.panel)

b0 <- 10
b1 <- 5
b2 <- .75
KX.raw <- K.raw%*%X
KX.spec <- Kspec%*%X
KX.raw.panel <- K.panel%*%X
KX.spec.panel <- Kspec.panel%*%X
WX.linraw <- Wlinraw%*%X
WX.invraw <- Winvraw%*%X
WX.linspec <- Wlinspec%*%X
WX.invspec <- Winvspec%*%X
WX.linraw.panel <- Wlinraw.panel%*%X
WX.invraw.panel <- Winvraw.panel%*%X
WX.linspec.panel <- Wlinspec.panel%*%X
WX.invspec.panel <- Winvspec.panel%*%X

###SLX w/ inflated DW
constant.DW <- c(-45, -25, 0, 25, 45)
par.est.SLX.DWcons.KX <- matrix(NA, nrow=reps, ncol=length(constant.DW)*6)
DW <- ifelse(D<50,50,D)
diag(DW) <- 0
DWX <- DW%*%X

for(K in 1:length(constant.DW)){
  DWhat <- DW + constant.DW[K]
  diag(DWhat) <- 0
  WX.hat <-DWhat%*%X
  for(i in 1:reps){
    Y <- b1*X + b2*DWX + rnorm(obs)
    model <- lm(Y~0 + X + WX.hat + KX.raw)
    vcv <- vcov(model)
    par.est.SLX.DWcons.KX[i,K*6-5] <- model$coef[1]
    par.est.SLX.DWcons.KX[i,K*6-4] <- model$coef[2]
    par.est.SLX.DWcons.KX[i,K*6-3] <- model$coef[3]
    par.est.SLX.DWcons.KX[i,K*6-2] <- sqrt(diag(vcv)[1])
    par.est.SLX.DWcons.KX[i,K*6-1] <- sqrt(diag(vcv)[2])
    par.est.SLX.DWcons.KX[i,K*6] <- sqrt(diag(vcv)[3])
  }
}


bias.SLX.DWcons.KX <- matrix(NA, nrow=4, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  bias.SLX.DWcons.KX[2,i] <- mean(par.est.SLX.DWcons.KX[,i*6-5]-b1)
  bias.SLX.DWcons.KX[3,i] <- mean(par.est.SLX.DWcons.KX[,i*6-4]-b2)
  bias.SLX.DWcons.KX[4,i] <- mean(par.est.SLX.DWcons.KX[,i*6-3])
}
bias.SLX.DWcons.KX

coverageprob.SLX.DWcons.KX <- matrix(NA, nrow=4, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  cp.beta1 <- coverage(par.est.SLX.DWcons.KX[,i*6-5], par.est.SLX.DWcons.KX[,i*6-2], b1, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX[2,i] <- cp.beta1$coverage.probability
  cp.b2 <- coverage(par.est.SLX.DWcons.KX[,i*6-4], par.est.SLX.DWcons.KX[,i*6-1], b2, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX[3,i] <- cp.b2$coverage.probability
  cp.Ky <- coverage(par.est.SLX.DWcons.KX[,i*6-3], par.est.SLX.DWcons.KX[,i*6], 0, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX[4,i] <- cp.Ky$coverage.probability
}
coverageprob.SLX.DWcons.KX

###SLX w/ inflated DW, panel
par.est.SLX.DWcons.KX.panel <- matrix(NA, nrow=reps, ncol=length(constant.DW)*8)
DW.panel <- ifelse(D.panel<50,50,D.panel)
DW.panel <- DW.panel*K.panel
DWX.panel <- DW.panel%*%X

for(K in 1:length(constant.DW)){
  DWhat <- DW.panel + constant.DW[K]
  DWhat <- DWhat*K.panel
  DWX.hat <-DWhat%*%X
  for(i in 1:reps){
    Y <- b0 + b1*X + b2*DWX.panel + rnorm(obs)
    model <- lm(Y~ X + DWX.hat + KX.raw.panel)
    vcv <- vcov(model)
    par.est.SLX.DWcons.KX.panel[i,K*8-7] <- model$coef[1]
    par.est.SLX.DWcons.KX.panel[i,K*8-6] <- model$coef[2]
    par.est.SLX.DWcons.KX.panel[i,K*8-5] <- model$coef[3]
    par.est.SLX.DWcons.KX.panel[i,K*8-4] <- model$coef[4]
    par.est.SLX.DWcons.KX.panel[i,K*8-3] <- sqrt(diag(vcv)[1])
    par.est.SLX.DWcons.KX.panel[i,K*8-2] <- sqrt(diag(vcv)[2])
    par.est.SLX.DWcons.KX.panel[i,K*8-1] <- sqrt(diag(vcv)[3])
    par.est.SLX.DWcons.KX.panel[i,K*8] <- sqrt(diag(vcv)[4])
  }
}

bias.SLX.DWcons.KX.panel <- matrix(NA, nrow=4, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  bias.SLX.DWcons.KX.panel[1,i] <- mean(par.est.SLX.DWcons.KX.panel[,i*8-7]-b0)
  bias.SLX.DWcons.KX.panel[2,i] <- mean(par.est.SLX.DWcons.KX.panel[,i*8-6]-b1)
  bias.SLX.DWcons.KX.panel[3,i] <- mean(par.est.SLX.DWcons.KX.panel[,i*8-5]-b2)
  bias.SLX.DWcons.KX.panel[4,i] <- mean(par.est.SLX.DWcons.KX.panel[,i*8-4])
}
bias.SLX.DWcons.KX.panel

coverageprob.SLX.DWcons.KX.panel <- matrix(NA, nrow=4, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  cp.beta0 <- coverage(par.est.SLX.DWcons.KX.panel[,i*8-7], par.est.SLX.DWcons.KX.panel[,i*8-3], b0, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX.panel[1,i] <- cp.beta0$coverage.probability
  cp.beta1 <- coverage(par.est.SLX.DWcons.KX.panel[,i*8-6], par.est.SLX.DWcons.KX.panel[,i*8-2], b1, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX.panel[2,i] <- cp.beta1$coverage.probability
  cp.rho <- coverage(par.est.SLX.DWcons.KX.panel[,i*8-5], par.est.SLX.DWcons.KX.panel[,i*8-1], b2, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX.panel[3,i] <- cp.rho$coverage.probability
  cp.KX <- coverage(par.est.SLX.DWcons.KX.panel[,i*8-4], par.est.SLX.DWcons.KX.panel[,i*8], 0, df=obs-model$rank)
  coverageprob.SLX.DWcons.KX.panel[4,i] <- cp.KX$coverage.probability
}
coverageprob.SLX.DWcons.KX.panel


###SAR with Inflated DW
par.est.SAR.DWcons.Ky.2SLS <- matrix(NA, nrow=reps, ncol=length(constant.DW)*8)
DWspec <- spectral(DW)

for(K in 1:length(constant.DW)){
  DWhat <- DW + constant.DW[K]
  diag(DWhat) <- 0
  DWhatspec <- spectral(DWhat)
  WX.hat <-DWhatspec%*%X
  for(i in 1:reps){
    y <- solve(I-b2*DWspec)%*%(X*b1 + rnorm(obs))
    Wy.DWhat <-DWhatspec%*%y
    Ky <- Kspec%*%y
    model <- tsls(y ~ 0+X + Wy.DWhat + Ky, ~0+X+WX.hat+KX.spec)
    model.lh <- summary(glht(model, linfct = c("Wy.DWhat + Ky = 0")))
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-7] <- model$coefficients[1]
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-6] <- model$coefficients[2]
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-5] <- model$coefficients[3]
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-4] <- model.lh$test$coefficients[1]
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-3] <- sqrt(model$V[1,1])
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-2] <- sqrt(model$V[2,2])
    par.est.SAR.DWcons.Ky.2SLS[i,K*8-1] <- sqrt(model$V[3,3])
    par.est.SAR.DWcons.Ky.2SLS[i,K*8] <- model.lh$test$sigma[1]
  }
}

bias.SAR.DWcons.Ky.2SLS <- matrix(NA, nrow=4, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  bias.SAR.DWcons.Ky.2SLS[1,i] <- mean(par.est.SAR.DWcons.Ky.2SLS[,i*8-7]-b1)
  bias.SAR.DWcons.Ky.2SLS[2,i] <- mean(par.est.SAR.DWcons.Ky.2SLS[,i*8-6]-b2)
  bias.SAR.DWcons.Ky.2SLS[3,i] <- mean(par.est.SAR.DWcons.Ky.2SLS[,i*8-5])
  bias.SAR.DWcons.Ky.2SLS[4,i] <- mean(par.est.SAR.DWcons.Ky.2SLS[,i*8-4]-b2)
}
bias.SAR.DWcons.Ky.2SLS

coverageprob.SAR.DWcons.Ky.2SLS <- matrix(NA, nrow=4, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  cp.beta1 <- coverage(par.est.SAR.DWcons.Ky.2SLS[,i*8-7], par.est.SAR.DWcons.Ky.2SLS[,i*8-3], b1, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS[1,i] <- cp.beta1$coverage.probability
  cp.b2 <- coverage(par.est.SAR.DWcons.Ky.2SLS[,i*8-6], par.est.SAR.DWcons.Ky.2SLS[,i*8-2], b2, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS[2,i] <- cp.b2$coverage.probability
  cp.Ky <- coverage(par.est.SAR.DWcons.Ky.2SLS[,i*8-5], par.est.SAR.DWcons.Ky.2SLS[,i*8-1], 0, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS[3,i] <- cp.Ky$coverage.probability
  cp.b2Ky <- coverage(par.est.SAR.DWcons.Ky.2SLS[,i*8-4], par.est.SAR.DWcons.Ky.2SLS[,i*8], b2, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS[4,i] <- cp.b2Ky$coverage.probability
}
coverageprob.SAR.DWcons.Ky.2SLS

###SAR w/ inflated DW, panel
par.est.SAR.DWcons.Ky.2SLS.panel <- matrix(NA, nrow=reps, ncol=length(constant.DW)*10)
DWspec.panel <- spectral(DW.panel)


for(K in 1:length(constant.DW)){
  DWhat <- DW.panel + constant.DW[K]
  DWhat <- DWhat*K.panel
  DWhatspec <- spectral(DWhat)
  DWX.hat <-DWhatspec%*%X
  for(i in 1:reps){
    y <- solve(I-b2*DWspec.panel)%*%(b0 + X*b1 + rnorm(obs))
    Wy.DWhat <-DWhatspec%*%y
    Ky <- Kspec.panel%*%y
    model <- tsls(y ~ X + Wy.DWhat + Ky, ~X+DWX.hat+KX.spec.panel)
    model.lh <- summary(glht(model, linfct = c("Wy.DWhat + Ky = 0")))
    vcv <- vcov(model)
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-9] <- model$coef[1]
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-8] <- model$coef[2]
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-7] <- model$coef[3]
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-6] <- model$coef[4]
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-5] <- model.lh$test$coefficients[1]
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-4] <- sqrt(diag(vcv)[1])
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-3] <- sqrt(diag(vcv)[2])
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-2] <- sqrt(diag(vcv)[3])
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10-1] <- sqrt(diag(vcv)[4])
    par.est.SAR.DWcons.Ky.2SLS.panel[i,K*10] <- model.lh$test$sigma[1]
  }
}


bias.SAR.DWcons.Ky.2SLS.panel <- matrix(NA, nrow=5, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  bias.SAR.DWcons.Ky.2SLS.panel[1,i] <- mean(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-9]-b0)
  bias.SAR.DWcons.Ky.2SLS.panel[2,i] <- mean(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-8]-b1)
  bias.SAR.DWcons.Ky.2SLS.panel[3,i] <- mean(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-7]-b2)
  bias.SAR.DWcons.Ky.2SLS.panel[4,i] <- mean(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-6])
  bias.SAR.DWcons.Ky.2SLS.panel[5,i] <- mean(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-5]-b2)
}
bias.SAR.DWcons.Ky.2SLS.panel

coverageprob.SAR.DWcons.Ky.2SLS.panel <- matrix(NA, nrow=5, ncol=length(constant.DW))
for (i in 1:length(constant.DW)){
  cp.beta0 <- coverage(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-9], par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-4], b0, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS.panel[1,i] <- cp.beta0$coverage.probability
  cp.beta1 <- coverage(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-8], par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-3], b1, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS.panel[2,i] <- cp.beta1$coverage.probability
  cp.rho <- coverage(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-7], par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-2], b2, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS.panel[3,i] <- cp.rho$coverage.probability
  cp.Ky <- coverage(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-6], par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-1], 0, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS.panel[4,i] <- cp.Ky$coverage.probability
  cp.rhoKy <- coverage(par.est.SAR.DWcons.Ky.2SLS.panel[,i*10-5], par.est.SAR.DWcons.Ky.2SLS.panel[,i*10], b2, df=obs-model$p)
  coverageprob.SAR.DWcons.Ky.2SLS.panel[5,i] <- cp.rhoKy$coverage.probability
}
coverageprob.SAR.DWcons.Ky.2SLS.panel

###Creating Tables

#Table 1
SLX.DWcons.KX.results <- matrix(NA,ncol=ncol(bias.SLX.DWcons.KX), nrow=2*nrow(bias.SLX.DWcons.KX))
colnames(SLX.DWcons.KX.results) <- c("-45", "-25", "0", "+25", "+45")
rownames(SLX.DWcons.KX.results)<- c("alpha", "","beta 1", "", "theta 1","","theta 2","")
for (i in 1:nrow(bias.SLX.DWcons.KX)){
  SLX.DWcons.KX.results[i*2-1,]<-bias.SLX.DWcons.KX[i,]
  SLX.DWcons.KX.results[i*2,]<-coverageprob.SLX.DWcons.KX[i,]
}

SLX.DWcons.KX.panel.results <- matrix(NA,ncol=ncol(bias.SLX.DWcons.KX.panel), nrow=2*nrow(bias.SLX.DWcons.KX.panel))
colnames(SLX.DWcons.KX.panel.results) <- c("-45", "-25", "0", "+25", "+45")
rownames(SLX.DWcons.KX.panel.results)<- c("alpha", "","beta 1", "", "theta 1","","theta 2","")
for (i in 1:nrow(bias.SLX.DWcons.KX.panel)){
  SLX.DWcons.KX.panel.results[i*2-1,]<-bias.SLX.DWcons.KX.panel[i,]
  SLX.DWcons.KX.panel.results[i*2,]<-coverageprob.SLX.DWcons.KX.panel[i,]
}

SLX.Results.2 <- rbind(SLX.DWcons.KX.results, SLX.DWcons.KX.panel.results)
SLX.Results.2 <- xtable(SLX.Results.2)
SLX.Results.2
print(SLX.Results.2, file="Table_1.txt")

#Table 2
SAR.DWcons.Ky.2SLS.results <- matrix(NA,ncol=ncol(bias.SAR.DWcons.Ky.2SLS), nrow=2*(nrow(bias.SAR.DWcons.Ky.2SLS)-1))
colnames(SAR.DWcons.Ky.2SLS.results) <- c("-45", "-25", "0", "+25", "+45")
rownames(SAR.DWcons.Ky.2SLS.results)<- c("beta 1", "", "Rho 1","","Rho 2","")
for (i in 1:(nrow(bias.SAR.DWcons.Ky.2SLS)-1)){
  SAR.DWcons.Ky.2SLS.results[i*2-1,]<-bias.SAR.DWcons.Ky.2SLS[i,]
  SAR.DWcons.Ky.2SLS.results[i*2,]<-coverageprob.SAR.DWcons.Ky.2SLS[i,]
}

SAR.DWcons.Ky.2SLS.panel.results <- matrix(NA,ncol=ncol(bias.SAR.DWcons.Ky.2SLS.panel), nrow=2*(nrow(bias.SAR.DWcons.Ky.2SLS.panel)-1))
colnames(SAR.DWcons.Ky.2SLS.panel.results) <- c("-45", "-25", "0", "+25", "+45")
rownames(SAR.DWcons.Ky.2SLS.panel.results)<- c("alpha", "", "beta 1", "", "Rho 1","","Rho 2","")
for (i in 1:(nrow(bias.SAR.DWcons.Ky.2SLS.panel)-1)){
  SAR.DWcons.Ky.2SLS.panel.results[i*2-1,]<-bias.SAR.DWcons.Ky.2SLS.panel[i,]
  SAR.DWcons.Ky.2SLS.panel.results[i*2,]<-coverageprob.SAR.DWcons.Ky.2SLS.panel[i,]
}

SAR.Results.2 <- rbind(SAR.DWcons.Ky.2SLS.results, SAR.DWcons.Ky.2SLS.panel.results)
SAR.Results.2 <- xtable(SAR.Results.2)
SAR.Results.2
print(SAR.Results.2, file="Table_2.txt")

SAR.Results.Lincom <-matrix(NA, ncol=ncol(bias.SAR.DWcons.Ky.2SLS.panel), nrow=4)
SAR.Results.Lincom[1,] <- bias.SAR.DWcons.Ky.2SLS[4,]
SAR.Results.Lincom[2,] <- coverageprob.SAR.DWcons.Ky.2SLS[4,]
SAR.Results.Lincom[3,] <- bias.SAR.DWcons.Ky.2SLS.panel[5,]
SAR.Results.Lincom[4,] <- coverageprob.SAR.DWcons.Ky.2SLS.panel[5,]
SAR.Results.Lincom


###Appendix
K <- matrix(1, nrow=obs, ncol=obs)
diag(K) <- 0
corr.Kx.Ux.results <- matrix(NA, nrow=reps, ncol=7)

for(j in 1:reps){
  a <- runif(1, 0, 1000)
  corr.Kx.Ux.results[j,1] <- a
  b <- runif(1, 0, 1000)
  corr.Kx.Ux.results[j,2] <- b
  c <- runif(1, -1000, 1000)
  corr.Kx.Ux.results[j,3] <- c
  x <- runif(obs, -a + c, a + c)
  Kx <- K%*%x
  U <- matrix(runif(obs*obs, -b, b), nrow=obs, ncol=obs)
  diag(U) <- 0
  Ux <- U%*%x
  corr <- cor.test(Kx, Ux)
  corr.Kx.Ux.results[j,4] <- corr$estimate
  corr.Kx.Ux.results[j,5] <- corr$conf.int[1]
  corr.Kx.Ux.results[j,6] <- corr$conf.int[2]
}

corr.Kx.Ux.results[,7] <- ifelse(corr.Kx.Ux.results[,5]<0 & corr.Kx.Ux.results[,6]>0, 0, 1)

corr.Kx.Ux.dataframe <- as.data.frame(corr.Kx.Ux.results)
summary(corr.Kx.Ux.dataframe$V4)
Average.Correlation <- density(corr.Kx.Ux.dataframe$V4) # returns the density data 
font_import()
loadfonts()
Cairo(file='Figure_1.pdf', type="pdf", dpi = 72)
plot(Average.Correlation, family="Verdana", main="Figure 1: Density Plot of the Correlation between Kx and Ux", xlab="")
dev.off() 
summary(corr.Kx.Ux.dataframe$V7) 



# Time Results
End.Time <- Sys.time()

End.Time - Start.Time
